library(readxl)
library(corrplot)
library(lm.beta)
library(Hmisc)
library(qgraph)
library(plyr)
library(dplyr)
library(psych)
library(neurobase)
library(oro.nifti)
library(dplyr)
library(png)
library(ggiraphExtra)
library(data.table)
library(cowplot)
library(ggplot2)
library(patchwork)
# library(NMF)
library(ggbeeswarm)
#source("/Volumes/Users/nfranzme/R_Projects/R_Imaging/workbench/workbench_map_volume_to_surface#.R")
library(ggpubr)
library(lubridate)
library(reticulate)
library(svMisc)
library(NetworkToolbox)
library(readxl)
library(writexl)
library(lme4)
library(lmerTest)
library(stringr)
library(boot)
library(table1)load input
define directories
dir.root = "/Volumes/NO NAME/fMRI_Day2_Graph_Theory/"
dir.fc.200 = paste0(dir.root, "/Schaefer200_functional_connectivity/")
dir.data.sheets = paste0(dir.root, "/data_sheets/")Data preparation
Step 1: load and prepare the data
The data frame contains data from the Alzheimer’s disease neuroimaging initiative (ADNI) database.
ADNI is a large multi-center study that focuses on multi-modal neuroimaging in Alzheimer’s disease patients, more information can be found here: http://adni.loni.usc.edu
Included imaging modalities in the current dataset are tau-PET, amyloid-PET and resting-state fMRI. We will subset the data to healthy controls (i.e. cognitively normal [CN] individuals without evidence for amyloid pathology [Ab.neg]), and subjects across the Alzheimer’s disease (AD) spectrum, defined as having abnormally elevated amyloid-beta (Ab) based on amyloid-PET. The AD spectrum includes individuals who are still cognitively normal (CN.Ab.pos, i.e. preclinical AD), subjects who show mild cognitive impairment (MCI.Ab.pos, i.e. prodromal AD), and patients with AD dementia (Dementia.Ab.pos)
For this course we will work primarily with cross-sectional data, stored in the “ADNI.data.bl” data frame. However, some individuals also have longitudinal tau-PET, you could have a look at those data to if there is time left.
-> read in data and select subjects across the AD spectrum
# read in ADNI data
ADNI.data <- readxl::read_xlsx(paste0(dir.data.sheets, "/ADNI_fMRI_PET.xlsx"))
ADNI.data <- subset(ADNI.data, ID %nin% c("sub-4431", "sub-2133", "sub-6559"))
# merge diagnosis and amyloid status
ADNI.data$DX.Ab <- paste0(ADNI.data$DX, ".", ADNI.data$amyloid.global.SUVR.status)
ADNI.data$DX.Ab <- factor(ADNI.data$DX.Ab,
levels = c("CN.Ab.neg",
"CN.Ab.pos",
"MCI.Ab.pos",
"Dementia.Ab.pos"))
# subset to controls and subjects across the AD spectrum (CN = cognitively normal, MCI = Mild Cognitive Impairment, Dementia = Alzheimer's disease dementia)
ADNI.data <- subset(ADNI.data, DX.Ab %in% c("CN.Ab.neg", "CN.Ab.pos", "MCI.Ab.pos", "Dementia.Ab.pos"))
# locate functional connectivity matrices
ADNI.data$FC.mat.200 <- str_replace(ADNI.data$FC.mat.200, "/Network/Cluster/ADNI/functional_connectivity/Schaefer_200/",
paste0(dir.root, "/Schaefer200_functional_connectivity/"))
# select baseline data
ADNI.data.bl = subset(ADNI.data, tau.pet.fu.yrs == 0)-> create some overview tables
First, you should get an idea of the type of data we’re working with. So we will check some basic distributions of biomarker and cognitive data for the current patient cohort.
create “table 1”
We can first create an overview table of demographics, as well as amyloid-PET and tau-PET load, to
see how biomarkers are distributed across diagnostic groups.
We’ll look at the following variables (Variable names on the left, explanation on the right)
age = age
sex = sex
centiloid = global amyloid-PET level (a typical “cut-off” for abnormality is 20)
tau.global.SUVR = global tau-PET level (a typical “cut-off” for abnormality is 1.3)
MMSE = Mini Mental State Exam, i.e. a global cognitive screening instrument, where lower values indicate stronger cognitive impairment
ADNI_MEM = A memory composite score across several memory tests. Lower values indicate stronger impairment
table1(~ age + factor(sex) + centiloid + tau.global.SUVR + MMSE + ADNI_MEM | DX.Ab, data = ADNI.data.bl)| CN.Ab.neg (N=192) |
CN.Ab.pos (N=67) |
MCI.Ab.pos (N=37) |
Dementia.Ab.pos (N=19) |
Overall (N=315) |
|
|---|---|---|---|---|---|
| age | |||||
| Mean (SD) | 72.2 (6.96) | 74.8 (5.96) | 75.9 (6.75) | 80.8 (8.09) | 73.7 (7.15) |
| Median [Min, Max] | 71.0 [56.0, 92.0] | 75.0 [65.0, 90.0] | 76.0 [61.0, 90.0] | 83.0 [62.0, 94.0] | 73.0 [56.0, 94.0] |
| factor(sex) | |||||
| F | 118 (61.5%) | 41 (61.2%) | 19 (51.4%) | 8 (42.1%) | 186 (59.0%) |
| M | 74 (38.5%) | 26 (38.8%) | 18 (48.6%) | 11 (57.9%) | 129 (41.0%) |
| centiloid | |||||
| Mean (SD) | -5.58 (12.8) | 65.6 (34.7) | 76.4 (31.7) | 93.0 (40.7) | 25.1 (45.6) |
| Median [Min, Max] | -7.12 [-34.3, 21.9] | 58.3 [23.2, 141] | 74.6 [22.8, 142] | 90.6 [44.2, 198] | 6.74 [-34.3, 198] |
| tau.global.SUVR | |||||
| Mean (SD) | 1.06 (0.0741) | 1.10 (0.0999) | 1.21 (0.222) | 1.36 (0.404) | 1.11 (0.164) |
| Median [Min, Max] | 1.06 [0.865, 1.45] | 1.10 [0.904, 1.62] | 1.15 [0.908, 1.83] | 1.23 [0.971, 2.37] | 1.08 [0.865, 2.37] |
| MMSE | |||||
| Mean (SD) | 29.3 (0.856) | 28.7 (1.74) | 26.7 (2.45) | 23.0 (4.73) | 27.7 (3.13) |
| Median [Min, Max] | 30.0 [27.0, 30.0] | 29.5 [24.0, 30.0] | 27.0 [22.0, 30.0] | 23.5 [14.0, 29.0] | 29.0 [14.0, 30.0] |
| Missing | 156 (81.3%) | 49 (73.1%) | 16 (43.2%) | 7 (36.8%) | 228 (72.4%) |
| ADNI_MEM | |||||
| Mean (SD) | 1.05 (0.566) | 0.924 (0.527) | 0.143 (0.584) | -0.685 (0.711) | 0.812 (0.740) |
| Median [Min, Max] | 0.991 [-0.229, 3.14] | 0.874 [-0.119, 1.94] | 0.0840 [-0.964, 1.69] | -0.509 [-2.39, 0.603] | 0.836 [-2.39, 3.14] |
| Missing | 1 (0.5%) | 0 (0%) | 0 (0%) | 0 (0%) | 1 (0.3%) |
As you can see, amyloid and tau levels increase (i.e. become more abnormal) across more severe
disease levels. Also, MMSE and ADNI_MEM values decrease, indicating stronger cognitive impairment.
plot biomarker distributions
Next, we should look at biomarker distributions of amyloid and tau pathology,
to get an idea of how amyloid and tau relate to increasing disease severity.
We’ll use ggplot2 for plotting. Some useful ggplot2 tutorials can be found here:
https://ggplot2.tidyverse.org
http://r-statistics.co/Complete-Ggplot2-Tutorial-Part1-With-R-Code.html
http://r-statistics.co/ggplot2-Tutorial-With-R.html
1. Amyloid-PET
Let’s first create a boxplot/beeswarm plot to illustrate the amyloid distribution across different diagnostic groups. The variables we’ll look at are centiloid (i.e. amyloid-PET) and DX.Ab (i.e. diagnostic groups)
# plot group differences
ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = centiloid)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 20, linetype = 2)
As you can see, amyloid levels increase across the AD spectrum. The dashed line indicates where amyloid becomes abnormal.
To statistically test these group differences, we can run an ANCOVA. It’s quite common standard to control for some usual covariates in clinical research, such as age and sex (and often also education)
So, lets quantify the group differences using an ANCOVA with a post-hoc Tukey Test. The ANCOVA tests whether there is an overall group-difference between diagnostic groups. The Post-Hoc Tukey test assesses the difference between single diagnostic groups.
# test group differences (ANCOVA), controlling for age and sex
## run ANCOVA
tmp.aov <- aov(data = ADNI.data.bl,
centiloid ~ DX.Ab + age + sex);
## summarize output
summary(tmp.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## DX.Ab 3 475451 158484 280.349 <2e-16 ***
## age 1 1435 1435 2.538 0.112
## sex 1 1078 1078 1.908 0.168
## Residuals 309 174681 565
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## run post-hoc Tukey test to determine group differences
TukeyHSD(tmp.aov, which = "DX.Ab")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = centiloid ~ DX.Ab + age + sex, data = ADNI.data.bl)
##
## $DX.Ab
## diff lwr upr p adj
## CN.Ab.pos-CN.Ab.neg 71.16379 62.4493273 79.87825 0.0000000
## MCI.Ab.pos-CN.Ab.neg 81.99277 70.9660878 93.01945 0.0000000
## Dementia.Ab.pos-CN.Ab.neg 98.54186 83.7714568 113.31227 0.0000000
## MCI.Ab.pos-CN.Ab.pos 10.82898 -1.7503318 23.40830 0.1191229
## Dementia.Ab.pos-CN.Ab.pos 27.37807 11.4151082 43.34104 0.0000770
## Dementia.Ab.pos-MCI.Ab.pos 16.54909 -0.7847542 33.88294 0.06734722. Tau-PET
Let’s do the same for tau-PET, so lets first create a boxplot/beeswarm plot to illustrate the tau distribution across diagnostic groups
# plot group differences
ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = tau.global.SUVR)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 1.3, linetype = 2)
As you can see, tau levels also increase across the AD spectrum. However, people usually start surpassing the threshold at symptomatic stages (MCI and Dementia). This is the case because amyloid precedes symptom onset by decades, while tau pathology is much closer to symptom onset. In fact, tau pathology is the key driver of neurodegeneration and cognitive decline, so you typically see abnormal tau first when people start to show symptoms.
To statistically test these group differences, we can again run an ANCOVA with a post-hoc Tukey Test
# test group differences (ANCOVA), controlling for age and sex
## run ANCOVA
tmp.aov <- aov(data = ADNI.data.bl,
tau.global.SUVR ~ DX.Ab + age + sex);
## summarize output
summary(tmp.aov)## Df Sum Sq Mean Sq F value Pr(>F)
## DX.Ab 3 1.998 0.6660 34.412 < 2e-16 ***
## age 1 0.350 0.3503 18.101 2.78e-05 ***
## sex 1 0.084 0.0841 4.345 0.0379 *
## Residuals 309 5.981 0.0194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## run post-hoc Tukey test to determine group differences
TukeyHSD(tmp.aov, which = "DX.Ab")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = tau.global.SUVR ~ DX.Ab + age + sex, data = ADNI.data.bl)
##
## $DX.Ab
## diff lwr upr p adj
## CN.Ab.pos-CN.Ab.neg 0.03695129 -0.01403989 0.08794248 0.2423433
## MCI.Ab.pos-CN.Ab.neg 0.15096161 0.08644085 0.21548237 0.0000000
## Dementia.Ab.pos-CN.Ab.neg 0.29530003 0.20887351 0.38172655 0.0000000
## MCI.Ab.pos-CN.Ab.pos 0.11401032 0.04040458 0.18761605 0.0004571
## Dementia.Ab.pos-CN.Ab.pos 0.25834874 0.16494414 0.35175334 0.0000000
## Dementia.Ab.pos-MCI.Ab.pos 0.14433842 0.04291236 0.24576449 0.00158463. Staging of tau pathology
One thing to keep in mind is that the above described analyses on tau-PET are based on “global” tau levels. However, tau does not - in contrast to amyloid - accumulate globally throughout the brain, but rather follows a consecutive spreading pattern that typically starts in the inferior temporal lobe. This “spreading pattern” of tau pathology has been already described in the nineties by Braak & Braak.
The surface rendering shows the six Braak-stages. Abnormal tau pathology is typically found first in Stage 1 (Entorhinal cortex) and then spreads to the hippocampus (Stage 2), the inferior temporal and limbic cortex (Stages 3&4) and finally to the association and primary somatosensory cortex (Stages 5&6). The plot on the left illustrates the likelihood of tau abnormality within each Braak stage across a group of AD patients, showing that early Braak stages typically show earlier tau abnormality than later Braak stages. Note that we usually discard the hippocampus when looking at tau-PET, because the tau-PET tracer shows quite significant off-target binding in the choroid plexus which is right next to the hippocampus. Thus, tau-PET signal in the hippocampus is quite confounded, at least for first-generation tau-PET tracers like flortaucipir.
Next, lets look at the distribution of tau pathology in the separate Braak stage ROIs. Tau-PET within these Braak-stage ROIs is labeled as tau.braak1.SUVR, tau.braak3.SUVR, tau.braak4.SUVR ,etc.
# plot group differences
# note that you can store plots in the workspace (e.g. p1 below) and later combine them using the patchwork package
p1 <- ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = tau.braak1.SUVR)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 1.3, linetype = 2) +
ggtitle("Braak 1")
p2 <- ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = tau.braak3.SUVR)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 1.3, linetype = 2) +
ggtitle("Braak 3")
p3 <- ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = tau.braak4.SUVR)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 1.3, linetype = 2) +
ggtitle("Braak 4")
p4 <- ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = tau.braak5.SUVR)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 1.3, linetype = 2) +
ggtitle("Braak 5")
p5 <- ggplot(data = ADNI.data.bl,
aes(
x = DX.Ab,
y = tau.braak6.SUVR)) +
geom_boxplot(alpha = 0.1, notch = T, width = 0.5) +
geom_beeswarm(alpha = 0.4) +
theme_minimal() +
geom_hline(yintercept = 1.3, linetype = 2) +
ggtitle("Braak 6")
# you can combine different plots using the patchwork package
# a "/" will add a new plot below, a "+" will add a plot next to the first one
p1 / p2 / p3 / p4 / p5
You can see that earlier Braak stage ROIs show abnormal signal more early during the disease course than later Braak stage ROIs
4. Association between amyloid and tau
Our current understanding of Alzheimer’s disease suggests that amyloid-beta deposition is the first pathological event to happen, that triggers downstream tau accumulation, neurodegeneration and cognitive decline. Thus, tau and amyloid should be correlated. Let’s have a look:
ggplot(ADNI.data.bl,
aes(
x = centiloid,
y = tau.global.SUVR
)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm") +
theme_minimal()## `geom_smooth()` using formula 'y ~ x'
This shows that higher amyloid is clearly associated with higher tau pathology. We can quantify the association using correlation or linear regression controlling for age and sex
# correlation
rcorr(ADNI.data.bl$centiloid,
ADNI.data.bl$tau.global.SUVR)## x y
## x 1.00 0.43
## y 0.43 1.00
##
## n= 315
##
##
## P
## x y
## x 0
## y 0# linear regression controlling for age and sex
tmp.lm <- lm(data = ADNI.data.bl,
tau.global.SUVR ~ centiloid + age + sex); summary(tmp.lm)##
## Call:
## lm(formula = tau.global.SUVR ~ centiloid + age + sex, data = ADNI.data.bl)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28412 -0.06856 -0.01542 0.04292 1.14154
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.3685946 0.0874300 15.654 < 2e-16 ***
## centiloid 0.0017479 0.0001875 9.320 < 2e-16 ***
## age -0.0040437 0.0012088 -3.345 0.000923 ***
## sexM -0.0198439 0.0168082 -1.181 0.238662
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.145 on 311 degrees of freedom
## Multiple R-squared: 0.2233, Adjusted R-squared: 0.2158
## F-statistic: 29.8 on 3 and 311 DF, p-value: < 2.2e-16TASK 1
Now that you know how to run group comparisons (i.e. ANCOVAs) and correlations/regression (linear models), try to check whether amyloid and tau are associated with cognition (e.g. MMSE and ADNI_MEM). How do amyloid and tau relate to cognitive performance? Try to create figures and statistics showing differences in cognitive performance between diagnostic groups, and associations between amyloid, tau and cognition.
fMRI analysis
Next lets go have a look at some fMRI parameters. We have prepared functional connectivity data using a common 200 ROI cortical atlas (see https://pubmed.ncbi.nlm.nih.gov/28981612/). A surface rendering of the atlas can be found below, including the seven networks (DMN = Default-Mode, DAN = Dorsal Attention, VAN = Ventral Attention, FPCN = Fronto-parietal control, Visual, Motor, Limbic)
1. Importing connectivity matrices
Atlas based connectivity is usually stored in matrix format also referred to as adjacency matrix. For 200 ROIs this means that the resulting adjacency matrix 200x200 elements large, where each element corresponds to the connectivity-metric between two ROIs (e.g. Fisher-z transformed Pearson correlation). All connectivity matrices are stored in the “Schaefer200_functional_connectivity” folder on the hard drive in .csv format. This folder is already stored in your R worksapce as “dir.fc.200”. The paths to the matrices are stored in the “ADNI.data.bl” data frame in the variable “FC.mat.200” Lets load an example to see how this looks like.
# import and format connectivity matrix using the read.csv command
current.adj.matrix = read.csv(ADNI.data.bl$FC.mat.200[1])When checking the dimension of this matrix, we realize that there are 201 columns and 200 rows. The first column contains the rownames, so we need to remove it in order to get a symmetrical 200x200 matrix
# check dimensions
dim(current.adj.matrix)## [1] 200 201# remove first column
current.adj.matrix.cleaned <- current.adj.matrix[,-1]
# clean diagnonal (i.e. remove infinite values with zeros)
diag(current.adj.matrix.cleaned) = 0
# transform to matrix format
current.adj.matrix.cleaned.matrix <- as.matrix(current.adj.matrix.cleaned)
# check dimensions
dim(current.adj.matrix.cleaned)## [1] 200 200Now the connectivity matrix has 200 rows and 200 columns. Lets have a look at the overall connectivity pattern. You can plot connectivity matrices using the “corrplot” package.
# plot correlation matrix
corrplot::corrplot(current.adj.matrix.cleaned.matrix,
diag = FALSE,
tl.pos = "n",
tl.cex = 0.5,
method = "color",
is.corr = FALSE)